# 2024/10/4

import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt

# Constants and parameters (fixed effects for PK and PD models)
Ka = 4.0489  # absorption rate constant
V = 2.561    # volume of distribution
Ke = 0.1439  # elimination rate constant
Tlag = 0.8454 # lag time for absorption

# PD parameters
Kgrow = 0.01    # tumor growth rate
Kkill = 0.1     # tumor killing rate
Kmax = 5        # maximum inhibition
KC50 = 300      # concentration for half-maximal killing rate
E0 = 350        # initial tumor volume

# Secondary parameter
TSC = Kgrow * KC50 / (Kkill * Kmax - Kgrow)

# Define the PK model
def pk_model(t, y, Ka, Ke):
    Aa, A1 = y
    dAa_dt = -Ka * Aa  # Absorption compartment equation
    dA1_dt = Ka * Aa - Ke * A1  # Central compartment equation
    return [dAa_dt, dA1_dt]

# Define the PD model
def pd_model(t, y, C, Kmax, KC50, Kkill, Kgrow):
    E1, E2, E3, E4 = y
    # Ensure non-negative tumor volumes
    E11 = max(E1, 0)
    E21 = max(E2, 0)
    E31 = max(E3, 0)
    E41 = max(E4, 0)
    
    # Inhibition function
    Inh = Kmax * (C / (KC50 + C))
    
    # Tumor growth and killing dynamics
    dE1_dt = (Kgrow - Inh * Kkill) * E11**(2/3)
    dE2_dt = (Inh * Kkill * E11**(2/3)) - Kkill * E21**(2/3)
    dE3_dt = Kkill * E21**(2/3) - Kkill * E31**(2/3)
    dE4_dt = Kkill * E31**(2/3) - Kkill * E41**(2/3)
    
    return [dE1_dt, dE2_dt, dE3_dt, dE4_dt]

# Initial conditions
AaDose = 1.0  # Initial dose
Aa0 = AaDose  # Initial amount of drug in absorption compartment
A10 = 0.1     # Initial amount of drug in central compartment
E1_0 = 0.7 * E0     # Initial tumor volume in E1
E2_0 = 0.2 * E0  # Initial tumor volume in E2
E3_0 = 0.1 * E0  # Initial tumor volume in E3
E4_0 = 0.0       # Initial tumor volume in E4

# Time points for simulation
t_span = (0, 100)  # Simulation from time 0 to 100 (arbitrary units)
t_eval = np.linspace(t_span[0], t_span[1], 1000)

# Solve the PK model
pk_sol = solve_ivp(pk_model, t_span, [Aa0, A10], args=(Ka, Ke), t_eval=t_eval)

# Drug concentration over time (C = A1 / V)
C_t = pk_sol.y[1] / V

# Solve the PD model using the concentration from PK model
def combined_model(t, y):
    # Interpolate C from the PK model solution at time t
    C = np.interp(t, pk_sol.t, C_t)
    return pd_model(t, y, C, Kmax, KC50, Kkill, Kgrow)

# Initial conditions for PD model
pd_sol = solve_ivp(combined_model, t_span, [E1_0, E2_0, E3_0, E4_0], t_eval=t_eval)

# Extract results
E1_t, E2_t, E3_t, E4_t = pd_sol.y
total_tumor_volume = E1_t + E2_t + E3_t + E4_t

# Plot results
plt.figure(figsize=(12, 6))

# Plot PK results
plt.subplot(1, 2, 1)
plt.plot(pk_sol.t, C_t, label="Concentration (C)")
plt.xlabel("Time")
plt.ylabel("Concentration (C)")
plt.title("PK Model")
plt.legend()

# Plot PD results
plt.subplot(1, 2, 2)
plt.plot(pd_sol.t, total_tumor_volume, label="Total Tumor Volume (E)")
plt.xlabel("Time")
plt.ylabel("Tumor Volume")
plt.title("PD Model")
plt.legend()

plt.tight_layout()
plt.show()